%%
tic

%%% This is the main file to solve the discrete time model.

clear all
close all
% load sv_discrete.mat
load initialV.mat

global parms i kr Hr jr cr Sf kf Nf cf gf sigmaf T0 lambdaf
y0 = 1.0;     % output of goods y0 - choose units so this is 1.0
%
Sbar = 2126.0527;    % Total feasible resources
Res_0 = 1./0.065100642;       % Initial estimate of proved reserves
rho_0 = Sbar - Res_0;    % Initial ratio of Alpha2 to Alpha3
g_0 = 0.11068716;      % Initial value of per unit mining cost
n0 = 0.00832913;      % Initial level of mining investment
%
% % gN_0 = -0.081639946;    % Partial derivative of g with respect to N at t=0
% % Alpha3 = -gS_0.*Rho_0./gN_0;
 
alpha3 = 15; %-gS_0.*Rho_0./gN_0;
alpha2 = rho_0.*alpha3;
%
delta = 0.04;
%
c0 = 0.6619974; % Initial value of consumption for calculate lambda_0
k0 = 3.6071282734; % Initial value of K for the differential equation
%
A =  y0./k0;
%
% Marginal cost of backstop energy p = (Gamma1+H)^(-alpha)
% Alpha is the slope of the learning curve
%
alpha = 0.25;
%
% Given alpha, Gamma1 determines the initial cost of renewable energy. Here we set it
% to 4 times the initial cost of fossil fuel.
%
Gamma1 = (4*g_0).^(-1./alpha);
%
% % After some time t, marginal cost will decline to Gamma2 and remain there.
% % We assume this ultimate minimum marginal cost of renewables
% % is 20% of the initial cost p when H = 0:
%
Gamma2 =0.8*g_0;
%
Abar = A*(1-Gamma2)+ (1-delta);
% Note: We should have Abar > 0. It is different from Abar in continuous
% time model

% To make discretized model comparable to continuous one, beta is redefined
% so that the long term growth rate of the two models are the same.
% (beta*Abar)^(1/gamma)-1 = 4.07%, from which beta = 0.9673
% beta_new = 1/(1+beta), where beta = 0.05.
beta  =  0.9524;

% psi is the effect of learning by doing on H
% psi has to be between 0 and 1
% A smaller value of psi will allow a larger role for explicit investment
% in renewable technology as opposed to learning by doing

psi = 0.33;
psr = 1/psi;
psc = 1-psi;
%
% Q = population growth that is used to "scale" resource extraction
% (variables other than fossil fuel exploitation are in per capita terms)
% popgr is the exogenous population growth rate

Q0 = 1;
popgr = 0.01;
%
% % gamma is the coefficient of relative rsik aversion. If we wish
% % to calibrate to a particular initial consumption level, we can
% % allow gamma to vary.
gamma = 4.0;
gamc = 1-gamma;
gamr = 1/gamma;

gS_0 = 0.00015;  % Partial derivative of g with respect to S at t=0
alpha1 = gS_0.*Res_0.^2;
alpha0 = g_0 - alpha1./Res_0;

% R&D and Learning spill-over. Expect that the former will be larger.
% thetaj = 0.5*psi;
% thetaB = 0.5*psi;
% thetaj = 0.05*psi;
% thetaB = 0.95*psi;
% thetaj = 0.95*psi;
% thetaB = 0.05*psi;
thetaj = 0;
thetaB = 0;

parms = [delta A Sbar alpha0 alpha1 alpha2 alpha3 Gamma1 alpha psi beta ...
    gamma Q0 popgr Abar];
% % Note that parms is a 1x14 vector



%%%%%%%%%%%%%%%  Regime 1: B>0, j>=0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% A very preliminary algorithm:
%%% 1. Given Tbar, decrease S_T0 first, until S_T0 is too small to solve the
%%% problem appropriately. Note that S_T0 has little effects on k0, but
%%% decreases N0 and S0, especially S0.
%%% 2. Increase k_Tbar until it is too large to solve the problem, which will
%%% increase k0, while decrease both N0 and S0.
%%% 3. If k0 is still smaller than the target, while N0 and S0 positive,
%%% decrease Tbar and repeat step 1 and 2 again.


%%% solution for guesses in continuous model
% Tbar = 315.8;
% k_Tbar = 141704.98998437249;
% S_T0 = 1613;
%%%
%1. Guess the time Tbar when the economy transits to the analytical
% regime. Tbar has no effects on renewable regime.

% 2. Guess the capital stock at k_Tbar.  On the
% other hand, the larger the k_Tbar is, the larger T0 is. k_T0 will be larger as well.

% 3. Guess S_T0, given Sbar = 2126.0527

% %%% PO case
Tbar = 366;
k_Tbar = 133937.035084899981;
S_T0 = 1595.42705305;
% T0 = 98, k0 = 3.607023,N0 = 9.2457e-05, S0 = 20.4151

% %%% rho = 0.5*psi, theta = 0.5*psi
% Tbar = 408;
% k_Tbar = 161864.00403179985;
% S_T0 = 1579.78478825;
% T0 = 107, k0 = 3.595431,N0 = -1.1463e-05, S0 = 20.3945


% %%% rho = 0.25*psi, theta = 0.75*psi
% k_Tbar is between 146197.46
% Tbar = 400;
% k_Tbar = 146197.463759450009;
% S_T0 = 1579.61845203;
% T0 = 106, k0 = 3.607139,N0 = 3.2377e-05, S0 = 20.4068


% %%% rho = 0.05*psi, theta = 0.95*psi
% Tbar = 387;
% k_Tbar = 136086.903981801035;
% S_T0 = 1582.48792067;
% T0 = 103, k0 = 3.608044,N0 = 8.9564e-05, S0 = 20.4153


% %%% rho(thetaj) = 0.95*psi, theta(B) = 0.05psi
% Tbar = 436; 
% k_Tbar = 204550.525554385968;
% S_T0 = 1570.19594116;
% T0 = 112, k0 = 3.601186,N0 = 9.5647e-05, S0 = 20.4126

% ======================================
% At Tbar, the marginal cost of backstop technology becomes Gamma2. Hence
% we have (Gamma1+H)^(-alpha)=Gamma2
H_Tbar = Gamma2^(-1/alpha)-Gamma1;

% Calculate k_t, H_t and V_t for all t in the renewable regime.

% Initialize the solution matrix
tt = 500; % tt should be larger than # of years in renewable regime
ir = zeros(tt,1);
jr = zeros(tt,1);
cr = zeros(tt,1);
kr = zeros(tt,1);
Hr = zeros(tt,1);
etar = zeros(tt,1);
lambdar = zeros(tt,1);
PRr = zeros(tt,1); %% p_t+1/p_t


while true
    % Value of variables at Tbar
    kr(1) = k_Tbar;
    Hr(1) = H_Tbar;
    ir(1) = ((beta*Abar)^gamr-1+delta)*k_Tbar;
    cr(1) = (1-Gamma2)*A*k_Tbar-ir(1);
    PRr(1) = 1/Abar;
    etar(1) = 0;
    lambdar(1) = cr(1)^(-gamma);
    i = 2;
    while true
        
        PRr(i) = 1/(A-(Gamma1+Hr(i-1))^(-alpha)*A+1-delta+(psi-thetaB)*jr(i-1)/((psc-thetaj)*kr(i-1)));
        cr(i) = (beta*cr(i-1)^(-gamma)/PRr(i))^(-gamr);
        lambdar(i) = cr(i)^(-gamma);
        Dbar = (psc-thetaj)*alpha*A*kr(i-1)*(Gamma1+Hr(i-1))^(-alpha-1)+A^(-psi)*kr(i-1)^(-psi)*jr(i-1)^psi;
        Cbar = Dbar*PRr(i);
        
        j = @(k) Cbar.^psr*A*k;
        inv = @(k) kr(i-1)-(1-delta)*k;
        H = @(k) Hr(i-1)-Cbar^(psr*psc)*A*k;
        
        if H(kr(i-1))>0 && inv(kr(i-1))>0
            kr(i) = fzero(@(k) (1-(Gamma1+H(k)).^(-alpha))*A.*k-inv(k)-j(k)-cr(i),kr(i-1));
            Hr(i) = H(kr(i));
            
            if Hr(i) == Hr(i-1)
                error('H converges to a positive number')
            end
            
            jr(i) = j(kr(i));
            ir(i) = inv(kr(i));
            etar(i) = beta*(etar(i-1) + lambdar(i-1)*alpha*(Gamma1+Hr(i-1))^(-alpha-1)*A*kr(i-1));
            i = i+1;
        else
            
            Hr(i) = 0;
            kr(i) = fzero(@(k) (1-Gamma1.^(-alpha))*A.*k-inv(k)-j(k)-cr(i),kr(i-1));
            jr(i) = j(kr(i));
            ir(i) = inv(kr(i));
            etar(i) = beta*(etar(i-1) + lambdar(i-1)*alpha*(Gamma1+Hr(i-1))^(-alpha-1)*A*kr(i-1));
            % Cut out the zero lines of solution matrix
            yearsr = i;
            Hr = Hr(1:yearsr);
            kr = kr(1:yearsr);
            jr = jr(1:yearsr);
            ir = ir(1:yearsr);
            cr = cr(1:yearsr);
            PRr = PRr(1:yearsr);
            etar = etar(1:yearsr);
            lambdar = lambdar(1:yearsr);
            break
        end
        
    end
    
    T0 = Tbar - yearsr;
    
    plot(T0+1:Tbar,jr(end:-1:1))
    xlabel('year')
    ylabel('Technology investment j')
    
    
    %%   %%%%%%%%%%%%%%%  Regime 2: R>0, n>=0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Key state variables: k and N
    % Control variables: i and n
    % Other state variables: Q and S
    % Transition point T0
    
    
    
    
    g = @(S,N) alpha0 + alpha1./(Sbar-alpha2./(alpha3+N)-S);
    gdS = @(S,N) alpha1*(alpha3+N).^2./((Sbar-S).*(alpha3+N)-alpha2).^2;
    gdN = @(S,N) -alpha1.*alpha2./((Sbar-S).*(alpha3+N)-alpha2).^2;
    
    while true
        
        T0 = Tbar - yearsr;
        Smax = Sbar;
        Smin = 0;
        
        % Set the initial results
        invf = zeros(T0+1,1);
        nf = zeros(T0+1,1);
        kf = zeros(T0+1,1);
        Nf = zeros(T0+1,1);
        Sf = zeros(T0+1,1);
        sigmaf = zeros(T0+1,1);
        cf = zeros(T0+1,1);
        Q = zeros(T0+1,1);
        gf = zeros(T0+1,1);
        PRf = zeros(T0+1,1);
        nuf = zeros(T0+1,1);
        lambdaf = zeros(T0+1,1);
        priceEf = zeros(T0+1,1);
        
        kf(1) = kr(end);
        cf(1) = cr(end);
        lambdaf(1) = lambdar(end);
        invf(1) = ir(end);
        PRf(1) = PRr(end);
        
        Q(1) = (1+popgr)^(T0);
        gf(1) = Gamma1^(-alpha)-psi*etar(end)*(A*kf(1))^(psi-1)*jr(end)^(1-psi)/(cf(1)^(-gamma));
        Sf(1) = S_T0;
        Nf(1) = alpha2/(Sbar-Sf(1)-alpha1/(gf(1)-alpha0))-alpha3;
        priceEf(1) = gf(1);
        
        if gf(1) < alpha0
            error(['g_T0 = ', num2str(gf(1)), ' g_T0 calculated from the renew regime is too small'])
            
        end
        
        z_T0 = [kf(1),Sf(1),Nf(1),lambdaf(1),sigmaf(1),nuf(1)];
        tspan0 = [T0, T0-1];
        % %
        % options = odeset('RelTol',1e-6,'AbsTol',1e-6,'events',@events0);
        % [tf,zf,TE0,ZE0,IE0] = ode45(@foss,tspan0,z_T0,options);
        %
        options = odeset('RelTol',5e-14,'AbsTol',5e-14);
        [tf,zf] = ode45(@foss,tspan0,z_T0,options);

        % Extract components of z for graphing
        
        kf(2)=zf(end,1);
        Sf(2)=zf(end,2);
        Nf(2)=zf(end,3);
        lambdaf(2)=zf(end,4);
        sigmaf(2)=zf(end,5);
        nuf(2)=zf(end,6);
        
        Q(2) = (1+popgr)^(T0-1);
        cf(2) = lambdaf(2)^(-gamr);
        gf(2) = g(Sf(2),Nf(2));
        invf(2) = inv(kf(2));
        PRf(2) = (beta*lambdaf(1))/lambdaf(2);
        priceEf(2) = gf(2)-sigmaf(2)*Q(2)/lambdaf(2);
        nf(2) = Nf(1)-Nf(2);
        
%             if priceEf(2)>priceEf(1)
%                warning(['abs(sigma*Q/lambda) too large: ', num2str(sigmaf(2)*Q(2)/lambdaf(2)) ])
%             end
        
        if lambdaf(2) - nuf(2) > 0
            disp(['T0 = ',num2str(T0), ', lambda = ',num2str(lambdaf(2)), ', nu = ', num2str(nuf(2))])
            %                 disp(['Decrease Tbar = ', num2str(Tbar), ' for 1'])
            %                 Tbar = Tbar - 1;
            
            disp(['Decrease S_T0 = ', num2str(S_T0), ' for 1'])
            Smax = S_T0;
            S_T0 = S_T0-1;
            continue
        end
        for i = 3:T0+1
            
            Q(i) = (1+popgr)^(T0-i+1);
            cf(i) = (beta*((1-gf(i-1)+sigmaf(i-1)*Q(i-1)/(cf(i-1)^(-gamma)))*A+1-delta))^(-gamr)*cf(i-1);
            lambdaf(i) = cf(i)^(-gamma);
            
            if ~isreal(cf(i))
                disp(['T0 = ', num2str(T0), ', cf(', num2str(i-1), ') = ', num2str(cf(i-1)), ', sigma*Q/lambda is too negative to make cf real'])
                %                         disp(['Increase Tbar = ', num2str(Tbar), ' for 1'])
                %                         Tbar = Tbar + 1;
                %
                disp(['Increase S_T0 = ', num2str(S_T0), ' for 1'])
                S_T0 = S_T0 + 1;
                %
                break
            end
            
            sigmaf(i) = beta*(sigmaf(i-1)-gdS(Sf(i-1),Nf(i-1))*A*kf(i-1)*lambdaf(i-1));
            
%             
%                     S =@(k) Sf(i-1)-Q(i)*A*k;
%                     N =@(k) cf(i)+kf(i-1)-(1-delta)*k+Nf(i-1)-A*k*(1-g(k,S(k)));
%                     kk = @(k) 1-g(S(k),N(k))+sigmaf(i)*Q(i)-delta/A+gdN(S(k),N(k))*k;
%                     kk0 = kf(i-1);
%                     options = optimset('Display','off','TolFun',1e-8,'TolX',1e-8);
%                     while true
%                         kf(i) = fzero(kk, kk0,options);
%                         if isnan(kf(i))
%                             if kk(kk0) < 0
%                             error(['When i = ', num2str(i), ', kk(kk0) = ', num2str(kk(kk0))])
%                             end
%                             while true
%                                 if kk0>10
%                                     kk0 = kk0-0.1;
%                                 elseif kk0 > 3
%                                     kk0 = kk0-0.01;
%                                 else
%                                     error('kf<3,too small')
%                                 end
%                                 if sign(kk(kk0)*kk(kf(i-1))) == -1;
%                                     break
%                                 end
%                                 
%                             end
%                         else
%                             break
%                         end
%                         
%                     end
%                     Nf(i) = N(kf(i));
            x0 = [kf(i-1),Nf(i-1)] ;
            options = optimset('Display','off','TolFun',1e-8,'TolX',1e-8);
            xopt = fsolve(@FOC_fossil,x0,options);
            kf(i) = xopt(1);
            Nf(i) = xopt(2);
            nf(i) = Nf(i-1)-Nf(i);
            if ~isreal(kf(i)) || isnan(kf(i)) || kf(i) < 0 || nf(i) <-1
                error(['Wrong kopt solution: ', 'kf(', num2str(i), ') = ', num2str(kf(i)), ', nf(', num2str(i), ') = ', num2str(nf(i))])
            end
            invf(i) = kf(i-1)-(1-delta)*kf(i);
            Sf(i) = Sf(i-1)-Q(i)*A*kf(i);
            gf(i) = g(Sf(i),Nf(i));
            % PRf(i) =1/(1-A*kf(i-1)*gdN(Sf(i-1),Nf(i-1)));
            PRf(i) = (beta*lambdaf(i-1))/lambdaf(i);
            priceEf(i) = gf(i)-sigmaf(i)*Q(i)/lambdaf(i);
            
            if i < T0 + 1 && Nf(i) < -1e-4
                disp(['When i=', num2str(i), ',N=',num2str(Nf(i)),', S=',num2str(Sf(i)), ', k=',num2str(kf(i))]);
                disp(['Increase S_T0 = ', num2str(S_T0), ' for 0.1']);
                if S_T0 <= Smax && S_T0 >= Smin
                    S_T0 = S_T0 + 0.1;
                else
                    disp(['Tbar is too large:', num2str(Tbar)]);
                    Smax = Sbar;
                    Smin = 0;
                    Tbar = Tbar-1;
                end
                break
                % error(['When i=', num2str(i), ',N=',num2str(Nf(i)),', S=',num2str(Sf(i)), ', k=',num2str(kf(i))]);
            end
            
            %     if Sf(i) < -1e-4
            %                 S_T0 = S_T0 + 0.1;
            %
            %                 error(['When i=', num2str(i), ',N=',num2str(Nf(i)),', S=',num2str(Sf(i)), ', k=',num2str(kf(i))]);
            %     end
        end
        %
        
        if i == T0+1
            disp(['N0 = ', num2str(Nf(end))]);
            if abs(Nf(end)) > 1e-4
                S_T0 = S_T0-0.1*Nf(end);
%                          elseif abs(Nf(end)) > 1e-3
%                              S_T0 = S_T0-0.01*Nf(end);
%                          elseif abs(Nf(end)) > 1e-4
%                              S_T0 = S_T0-0.01*Nf(end);
            else
                break
            end
        end
        
    end
    disp(['Tbar = ',num2str(Tbar),';']);
    disp(['k_Tbar = ', num2str(k_Tbar,18), ';']);
    disp(['S_T0 = ',num2str(S_T0,12),';']);
    disp(['T0 = ', num2str(T0),', k0 = ', num2str(kf(end),7), ',N0 = ',num2str(Nf(end)),', S0 = ',num2str(Sf(end))]);
    %
    kdiff = kf(end)-k0;
    disp(['kdiff = ', num2str(kdiff,8)]);
    if abs(kdiff) > 0.0001
        k_Tbar = k_Tbar + 1e-5*kdiff;
    else
        break
    end
%   break
end

PR = [PRf(end:-1:1);PRr(end:-1:1)];
price = ones(Tbar+1,1);

for i = 2: Tbar+1
    price(i) = price(i-1)*PR(i);
end

% Compare the results of the two models.

figure(1)
subplot(2,2,1), plot(T0+1:Tbar,kr(end:-1:1))
xlabel('year')
ylabel('')
title('(a) Capital stock k')

subplot(2,2,2),  plot(T0+1:Tbar,Hr(end:-1:1))
xlabel('year')
ylabel('')
title('(b) Cumulative knowledge H')

subplot(2,2,3),  plot(T0+1:Tbar,ir(end:-1:1))
xlabel('year')
ylabel('')
title('(c) Capital investment i')

subplot(2,2,4),  plot(T0+1:Tbar,jr(end:-1:1))
xlabel('year')
ylabel('')
title('(d) Technology investment j')

figure(2)
subplot(2,3,1), plot(kf(end:-1:1))
xlabel('year')
ylabel('')
title('(a) Capital stock k')
subplot(2,3,2), plot(Nf(end:-1:1))
xlabel('year')
ylabel('')
title('(b) Techonology progress N')
subplot(2,3,3), plot(Sf(end:-1:1))
xlabel('year')
ylabel('')
title('(c) Cumulative fuel use S')
subplot(2,3,4), plot(invf(end:-1:1))
xlabel('year')
ylabel('')
title('(d) Capital investment i')
subplot(2,3,5), plot(nf(end:-1:1))
title('(e) Technology investment n')
xlabel('year')
ylabel('')
subplot(2,3,6), plot(gf(end:-1:1))
xlabel('year')
ylabel('')
title('(f) Mining cost g')

figure(3)
subplot(2,3,1), plot(sigmaf(end:-1:1))
xlabel('year')
ylabel('')
title('sigma')
subplot(2,3,2), plot(priceEf(end:-1:1))
xlabel('year')
ylabel('')
title('energy price')
subplot(2,3,3), plot(lambdaf(end:-1:1))
xlabel('year')
ylabel('')
title('lambda')
subplot(2,3,4), plot(lambdaf(end:-1:1).*priceEf(end:-1:1))
xlabel('year')
ylabel('')
title('nominal energy price')
subplot(2,3,5), plot(Q(end:-1:1))
title('population')
xlabel('year')
ylabel('')
subplot(2,3,6), plot(sigmaf(end:-1:1)/lambdaf(end:-1:1))
xlabel('year')
ylabel('')
title('sigma/lambda')

grlr = (beta*Abar)^(gamr)-1;
gr_cf = (cf(1)/cf(end))^(1/T0)-1;
gr_kf = (kf(1)/kf(end))^(1/T0)-1;
gr_cr = (cr(1)/cr(end))^(1/(Tbar-T0))-1;
gr_kr = (kr(1)/kr(end))^(1/(Tbar-T0))-1;
toc
